NEWTON'S METHOD AND SYMMETRY FOR SEMILINEAR ELLIPTIC PDE 

ON THE CUBE 



JOHN M. NEUBERGER, NANDOR SIEBEN, AND JAMES W. SWIFT 



Abstract. We seek discrete approximations to solutions u : Q. — s> M of semilinear elliptic partial 
differential equations of the form Alt + f s (u) = 0, where f s is a one-parameter family of nonlinear 
C*") ■ functions and Q is a domain in R d . The main achievement of this paper is the approximation of 

solutions to the PDE on the cube Q = (0, 7r) 3 C R 3 . There are 323 possible isotropy subgroups 
of functions on the cube, which fall into 99 conjugacy classes. The bifurcations with symmetry 
in this problem are quite interesting, including many with 3-dimensional critical eigenspaces. Our 
automated symmetry analysis is necessary with so many isotropy subgroups and bifurcations among 
them, and it allows our code to follow one branch in each equivalence class that is created at a 
bifurcation point. Our most complicated result is the complete analysis of a degenerate bifurcation 
Q\ , with a 6-dimensional critical eigenspace. 

This article extends the authors' work in Automated Bifurcation Analysis for Nonlinear Elliptic 
Partial Difference Equations on Graphs (Int. J. of Bifurcation and Chaos, 2009), wherein they 
combined symmetry analysis with modified implementations of the gradient Newton-Galerkin algo- 
rithm (GNGA, Neuberger and Swift) to automatically generate bifurcation diagrams and solution 
graphics for small, discrete problems with large symmetry groups. The code described in the cur- 
rent paper is efficiently implemented in parallel, allowing us to investigate a relatively fine-mesh 
discretization of the cube. We use the methodology and corresponding library presented in our 
paper An MPI Implementation of a Self-Submitting Parallel Job Queue (Int. J. of Parallel Prog., 
2012). 



C/3 
Qh 



> 

in 

00 

o 
o 



1. Introduction 

We are interested in finding and approximating solutions u : £7 — > R of semilinear elliptic equa- 
tions with zero Dirichlet boundary conditions, 



(TV J Au + f s (u) = in O 



u = on dil, 

where f s : R — > R satisfies / s (0) = 0, f' s (0) = s, and Q is a region in M? or M 3 . Our code also 
works for zero Neumann boundary conditions, and a wide range of nonlinearities. In this paper we 
present results for PDE (pQ) with f s (t) = st + i 3 on the square SI = (0, ir) 2 , and more challengingly, 
on the cube SI = (0, 7r) 3 . By finding and following new, bifurcating branches of (generally) lesser 
symmetry we are able to approximate, within reason, any solution that is connected by branches to 
the trivial branch. The more complicated solutions bifurcating farther from the origin (u, s) = (0, 0) 
are of course progressively more challenging to locate and accurately approximate. 

Generally, we apply Newton's method to the gradient of an action functional whose critical 
points are solutions to our PDE. For an exposition of our initial development of the gradient 
Newton-Galerkin algorithm (GNGA) and our first application of it to the square, see |16| . 

This article extends the methods for so-called partial difference equations (PdE) from [13] to 
large graphs, that is, fine mesh discretizations for PDE. For small graphs with possibly large 
symmetry groups, the code in [13] automated the analysis of symmetry, isotypic decomposition, 
and bifurcation. We use here the GAP (Groups, Algorithms, and Programming, see [8]) and 
Mathematica codes presented in those articles to automatically generate a wealth of symmetry 
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information for use by our branch-following C++ code. Some of this information is summarized 
in the bifurcation digraph, which shows the generic symmetry-breaking bifurcations. In |13] we 
developed two modified implementations of the gradient Newton-Galerkin algorithm, namely the 
tangent algorithm (tGNGA) for following bifurcation curves and the cylinder algorithm (cGNGA) 
for switching branches at bifurcation points. Together with a secant method for locating these 
bifurcation points, in the current PDE setting we are able to handle most difficulties that arise 
when encountering accidental degeneracies and high-dimensional critical eigenspaces. 

Since we use here a fine mesh to investigate PDE (JTJ) on the cube O = (0,7r) 3 , the practical 
implementation of the algorithms from [13] requires increased efficiencies. In particular, we use 
isotypic decompositions of invariant fixed-point spaces to take advantage of the block diagonal 
structure of the Hessian matrix. In doing so, we substantially reduce the dimension of the Newton 
search direction linear system, and hence also reduce the number of costly numerical integrations. 
The same theory allows for reduced dimensions in many of the search spaces when seeking new, 
bifurcating solutions near high-dimensional bifurcation points in the presence of symmetry. 

Even with these efficiency improvements, we found it necessary to convert our high-level branch 
following strategy to use parallel computing. Most of the details of the parallel implementation 
can be found in [15], where we present a general methodology using self-submitting job queues to 
implement many types of mathematical algorithms in parallel. In particular, therein we develop a 
light-weight, easy-to-use C++ parallel job queue library, which we have used in obtaining the cube 
results found in this article. 

Our numerical results are summarized in bifurcation diagrams, which plot a scalar function such 
as the value of u at a generic point u(x* ,y* , z*) versus s, for approximate solutions to Equation ([1]) 
with parameter s. These diagrams indicate by line type the Morse Index (MI) of solutions, which 
typically changes at bifurcation and turning points. We present graphics for individual approximate 
solutions in several formats. For the most part, we find that representative graphics using a small, 
fixed collection of patches ("flags") most clearly show the symmetries of real- valued functions 
of three variables. We call these flag diagrams. We also include some contour plots of actual 
solution approximations. A fairly comprehensive collection of graphics and supporting information 
describing the symmetries of functions on the cube, possible bifurcations of nonlinear PDE on the 
cube, and more example approximate PDE solutions can be found on the companion website [14] . 

In [13], we considered many small graphs where scaling was not used, and hence PDE were 
not involved. In the present setting, we approximate a solution u to PDE (pQ) with a vector 
u = (u„) € M. N whose components represent u values at N regularly spaced grid points in Q 
located a distance Ax apart. Thus, our approach is equivalent to applying our algorithms to the 
finite dimensional semilinear elliptic partial difference equation (PdE) 

(2) -Lu + / s (u) = mQ N , 

where O^y is a graph with N vertices coming from a grid. The matrix L is in fact the graph 
Laplacian on Qjy, scaled by /^ x y2 , and modified at boundary vertices to enforce a zero-Dirichlet 
problem boundary condition. See |11] for a discussion of ghost points for enforcing boundary 
conditions. For general regions in R 2 and M 3 , we approximate eigenvectors of L using standard 
linear techniques, e.g., Matlab's eigs or some other easy to use implementation of ARPACK. For the 
square and cube the eigenfunctions are of course well known explicitly in terms of sine functions, 
so the consideration of L is not necessary. Since accurate PDE results require the dimension N to 
be very large, in the expansions of our approximate solutions u = Ylm=i a rnM> m we use M <^ N 
discretized eigenfunctions of —A with this boundary condition. 

In Section[2]we present some theory for the action functional, its gradient and Hessian, symmetry 
of functions, the corresponding fixed-point subspaces, and bifurcations with symmetry. We apply 
the general symmetry theory to the basis generation process for the square and cube. In Section [3] 
we outline the algorithms used in our project. We include a high-level description of our numerical 
methods and corresponding implementations, including the new use of self-submitting parallel job 
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queues applied to obtain accurate high-resolution solutions for the cube. We also describe our 
method for taking advantage of the block structure of the Hessian and our procedure for generating 
contour plots of approximate solutions. Section [4] contains the results from our experiment on the 
square, essentially an efficient and automatic refinement of the computations found in |16j . 

Our main numerical results are found in Section [5j Namely, we investigate PDE (TjQ) on the 
cube, where it is required to use a large number of eigenmodes and spacial grid points in order 
to find nodally complicated solutions of high MI, lying in many different fixed-point spaces of the 
fairly large symmetry group. We present several interesting examples from the companion website. 
The website shows examples of a solution with each of the symmetry types that we found in our 
investigation. Due to space limitations, we do not show all of these solutions in this paper. Rather, 
we concentrate on the first six primary bifurcations, and one of the secondary bifurcations with 
symmetry (the symmetry of a tetrahedron) . Two of the primary bifurcations that we consider have 
degenerate bifurcation points, including one with a six-dimensional eigenspace that is the direct 
sum of two 3-dimensional irreducible representations of Oh, the symmetry of the cube. 

There does not seem to be much in the literature that specifically investigates the bifurcation 
and symmetry of solutions to semilinear elliptic PDE on the cube. The article [3] is interesting for 
pushing the nonlinearity power to the critical exponent in the cube case. The interested reader 
can consult works by Zhou and co-authors for alternate but related methods and algorithms for 
computing solutions to semilinear elliptic PDE, e.g., [18 } 1191 120] and the recent book [5]. 

2. Symmetry and Invariance 

For the convenience of the reader, in this section we summarize enough notation and theory 
from [13j to follow our new results. We also include square and cube-specific information required 
to apply our algorithms in our particular cases. 

2.1. The Functional Setting. Our techniques rely on two levels of approximation, namely the 
restriction of functions to a suitably large M-dimensional Galerkin subspace Bm Q H = Hq(JI), 
and the discretization of 0, to £1^. We call the natural numbers M and N the Galerkin and spacial 
dimensions of our approximations, respectively. For the regions considered in this paper, it suffices 
to divide the region up into N squares or cubes with edge length Ax, and then place a gridpoint 
x n in the center of each such cell. The graph Ojy has a vertex Vi corresponding to each gridpoint, 
with edges e n j determined by the several neighbors Xj which are at distance Ax away from x n . 
With this arrangement, the simple numerical integration scheme used below in Equations [6] and [8] 
to evaluate the nonlinear terms in our gradient and Hessian computations becomes the midpoint 
method. 

The eigenvalues of the negative Laplacian with the zero Dirichlet boundary condition satisfy 
(3) < Ai < A 2 < A 3 < ► oo, 

and the corresponding eigenfunctions {V>m}meN can be chosen to be an orthogonal basis for the 
Sobolev space H, and an orthonormal basis for the larger Hilbert space 1? = L 2 (Q), with the usual 
inner products. In the cases Q = (0, 7r) d for d = 2 and d = 3, we take the first M eigenvalues 
(counting multiplicity) from 

{Xij := i 2 +j 2 \i,j <E N} or {X ijj)k := i 2 + j 2 + k 2 | k £ N}, 

and singly index them in a vector A = (Ai, . . . , Xm)- I n these cases, the corresponding eigen- 
functions Vm that we use are appropriate linear combinations of the well-known eigenfunctions 

i^i,j(x,y) = f sin(ix) sin(jy) and ipij t k(x, y, z) = (f) 3//2 sin(ix) sin(jy) sm(kz). We process the 
eigenfunctions using the projections given in Section 12.31 in order to understand and exploit the 
symmetry of functions in terms of the nonzero coefficients of their eigenfunction expansions. The 
ip m are discretized as tp m € K , m € {1, . . . ,M}, by evaluating the functions at the gridpoints, 
i.e., i\) m = (tp m (xi), . . . ,ip m (x]y)). The eigenvectors ip m form an orthonormal basis for an M- 
dimensional subspace of R . 
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Using variational theory, we define a nonlinear functional J whose critical points are the solutions 
of PDE ([T|). We use the Gradient-Newton-Galerkin- Algorithm (GNGA, see [16 i) to approximate 
these critical points, that is, we seek approximate solutions u lying in the subspace 

B M ■■= span{V>i, . . . , if) M } C H, 

which in turn are discretely approximated in R^ by 

M 

m=l 

The coefficient vectors a in R (simultaneously the approximation vectors u in R ) are computed 
by applying Newton's method to the eigenvector expansion coefficients of the approximation —L + 
/ s (u) of the gradient VJ s (u). 

Let F s (p) = Jq f s (t) dt for all p G R define the primitive of f s . We then define the action 
functional J s : H — > R by 

M 

(4) J s (u) 



f i | V n| 2 - F s {u) dV=J2 \a 2 m \ m - [ F s (u) dV. 
Jn m=1 Jn 



The class of nonlinearities f s found in [H Uj imply that J s is well defined and of class C 2 on H. 
Computing directional derivatives and integrating by parts gives 

(5) J' s (u)(ip m ) = / Vu ■ VV> m - f s (u) ipm dV = a m \ m - / f s (u) ipm dV , 

Jn Jn 

for m £ {1, . . . , M}. Replacing the nonlinear integral term with a sum that is in fact the midpoint 
method given our specific (square or cube) grid gives the gradient coefficient vector g G R M defined 
by 

N 

(6) g m = a m X m - ^2 fs(Un)(ll> m )n AV. 

n=l 

Here, the constant mesh area or volume factor is given by AV = Vol(Q)/N = n d /N. The functions 
Pb m VJ s (u) and Ylm=i9m''Pm are approximately equal and are pointwise approximated by the 
vector — Lu + / s (u). 

To apply Newton's method to find a zero of g as a function of a, we compute the coefficient 
matrix h for the Hessian as well. A calculation shows that 

(7) J'^u)(M m ) = fvin- Vipm - f' s (u) rPi ip m dV = \ t 5 lm - [ f s (u) ^ ^ m dV, 

Jn Jn 

where 5i m is the Kronecker delta function. Again using numerical integration, for I, m G {1, . . . , M} 
we compute elements of h by 

N 

(8) hi m = Xidim ~ ^2 fs{u n )('4>l)n{'4 } m)n AV. 

n=l 

The coefficient vector g G R and the M x M coefficient matrix h represent suitable projections 
of the L 2 gradient and Hessian of J, restricted to the subspace Bm, where all such quantities are 
defined. The least squares solution \ to the M-dimensional linear system h\ = g always exists 
and is identified with the projection of the search direction (D\ J s (u)) _1 V2«/ S (w) onto Bm- The L2 
search direction is not only defined for all points u G Bm such that the Hessian is invertible, but is 
in that case equal to (D 2 H J s (n)) _1 V// J s (u). 

The Hessian function h s : R A/ R M or h : R M+1 -)■ R M is very important for identifying 
bifurcation points. If h(a* , s*) is invertible at a solution (a*, s*), then the Implicit Function Theorem 
guarantees that there is locally a unique solution branch through (a*,s*). When h(a*,s*) is not 
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Figure 1. Small graphs used to generate symmetry information for analyzing 
functions on the square and cube, respectively. The graph on the right, with full 
octahedral symmetry, is the skeleton graph of the great rhombicuboctahedron, which 
has 48 vertices, 72 edges, and 26 faces. This solid can be inscribed in the cube, 
with 8 vertices in each of the 6 faces. Thus, there is a one-to-one correspondence 
between elements of and vertices of the skeleton graph, after one vertex (chosen 
arbitrarily) has been assigned to the identity element. 

invertible, then (a*,s*) is a candidate for a bifurcation point, defined in the next subsection. The 
kernel 

E = Null h(a*,s*) 

of the Hessian at a bifurcation point is called the critical eigenspace. A Lyapunov-Schmidt reduction 
of the gradient g{a* , s*) can be done to obtain the Lyapunov-Schmidt reduced gradient g : E — > E 
[9]. Local to the point where h is singular, there is a one-to-one correspondence between zeros 
of g and zeros of g. The Lyapunov-Schmidt reduced bifurcation equations are g = 0. We refer to 
the reduced gradient or reduced bifurcation equations when the Lyapunov-Schmidt reduction is 
understood. 

Newton's method in coefficient space is implemented by fixing s, initializing the coefficient vector 
a with a guess, and iterating 

(9) a ^— a — x, where h(a, s)x = g(a, s). 

When it converges, the algorithm converges to vectors a and u = J2m=i a m^ m giving g = 0, and 
hence an approximate solution to PDE ([I]) has been found. The search direction \ in Newton's 
method is found by solving the system in ([9]) without inverting h(a, s). The solver we use returns 
the least squares solution for an overdetermined system and the solution with smallest norm for an 
underdetermined system. We observe experimentally that Newton algorithms work well even near 
bifurcation points where the Hessian is not invertible. In Section [3] we include brief descriptions 
of the tGNGA and cGNGA, the modifications of the GNGA actually implemented in our current 
code. 

In Sections 12.31 and 12.41 we explain the details of our method for constructing a specific basis 
of eigenfunctions that allows our code to exploit and report symmetry information. Regardless 
of whether we know a basis for Bm in terms of sines or must approximate one using numerically 
computed eigenvectors of the sparse matrix L, we are required to first compute various symmetry 
quantities relevant to the region 0, and the possible symmetry types of expected solutions. For 
this we use GAP [8]. We start with a graph G which has the same symmetry as fijv> but with 
a significantly smaller number of vertices. In Figure [1] we show a 12-vertex graph with the E84 
symmetry of the square, and a 48-vertex graph with the Q/j symmetry of the cube. These graphs 
are the smallest we found that allow functions on the vertices to have all possible symmetry types. 
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For example, the square with four vertices does not allow a function with only four-fold rotational 
symmetry. The group information used to analyze the symmetry of all functions on the square 
and cube regions was done automatically by our suite of programs. The automated GAP processes 
adapted from [13] would fail if fed instead the exceedingly large graphs Qn- Our GAP code is 
applied to the smaller graph G to generate all the files encoding the bifurcation digraph and the 
underlying fixed-point space decompositions. The results from the small-graph analysis inform us 
and our code about the symmetries of functions on f2 and the vectors in M. N which approximate 
such functions at the gridpoints. 

2.2. Symmetry of functions. We assume that f s is odd; the case where f s is not odd is easily 
inferred. To discuss the symmetry of solutions to Equation (pQ), we note that Aut(f2) x Z 2 = 
Aut(f27\r) x Z2 = Aut(G) x Z2, where G is one of the small graphs depicted in Figured! We define 

r = Aut(^Ar) x Z 2 , 

where Z2 = {1,-1} is written multiplicatively. 
The natural action of T on R N is defined by 

(10) (7-u)i = /3w 7r -i W , 

where 7 = (a w ,/3) G Fq and u G R . We usually write a for (a, 1) and —a for (a, — 1). The 
symmetry of u is the isotropy subgroup Sym(u) := Stab(u, Tq) = {7 G Tq | 7 • u = u}. The 
symmetries Sym(u) of functions u : f2 — > R are isomorphic to Sym(u). Two subgroups Tj and Tj of 
To are called conjugate if Tj = 7rj7 _1 for some 7 G Tq. The symmetry type of u is the conjugacy 
class [Sym(u)] of the symmetry of u; a similar definition holds for the symmetry type of a function 
u : 0, — > R. We say that two symmetry types are isomorphic if they have isomorphic representatives. 
We use the notation Q := {Tq, . . . , T q } for the set of symmetries and S := {Sq = [Tq], . . . , S r } for 
the set of symmetry types. 

Let X be the set of all solutions (u, s) to (|2j) in M, N x R. We define a branch of solutions to 
be a maximal subset of X that is a C 1 manifold with constant symmetry. The trivial branch 
{(0, s) I s G M} contains the trivial solution u = 0, which has symmetry Tq if f s is odd, and 
symmetry Aut(G) otherwise. A bifurcation point is a solution in the closure of at least two different 
solution branches. We call the branch containing the bifurcation point the mother, and the other 
branches, for which the bifurcation point is a limit point, are called daughters. Note that there is 
not a bifurcation at a fold point, where a branch of constant symmetry is not monotonic in s. 

The action of Tq on 1^ induces an action of Tq on R , given the correspondence of functions 
and coefficient vectors u = Ylm=i V>m- The gradient function g s : M A/ — > R M is To-equivariant, i.e., 
g s (j-a) = j-g s (a) for all 7 G Tq, a G M M , and s G R. As a consequence, if (a, s) is a solution to ([2]) 
then (7 • a, s) is also a solution to ([2]), for all 7 G Tq. Following the standard treatment [13], for 
each Ti < Tq we define the fixed-point subspace of the Tq action on M M to be 

Fix(r i; R M ) = {a G R M | 7 ■ a = a for all 7 G Ti}. 

The fixed-point subspaces for any of the function spaces V = Gm, R n > E, or H is defined as 

Tix(Ti, V) = {u G V I 7 • a = a for all 7 G Ti}. 

There is a one-to-one correspondence between o G M M and n G Ga/> and we will often write Fix(Fj) 
when the equation is valid for any ambient space. These fixed-point subspaces are important 
because they are g s -invariant, meaning that g s (Fix(Fj)) C Fix(Fj). For efficiency in our code, we 
restrict g s to one of these fixed-point subspaces, as described in Subsection 13.21 

As mentioned in the previous subsection, the reduced bifurcation equations are g = 0, where 
g : E — > E is the reduced gradient on the critical eigenspace E. A fold point (a*,s*) is not a 
bifurcation point even though h(a*,s*) is singular. 

Consider a bifurcation point (a*,s*) in coefficient space, or the corresponding (u*,s*), and let 
Ti = Sym(u*) be the symmetry of the mother solution. Then there is a natural action of Ti on 
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E, and the reduced gradient g is equivariant. That is, 5(7 • e,s) = 7 • g(e,s) for all 7 G Tj, and 
e £ E. Let T- be the kernel of the action of Ti on E. Then T^/T^ acts freely on and we say that 
the mother branch undergoes a bifurcation with T^/T^ symmetry. This bifurcation is generic, or 
non- degenerate, if the action of Tj/r^ on E is irreducible, and other non-degeneracy conditions (see 
[9]) are met. For each Tj < Ti, if e 6 Fix(r.,-, £7), then Sym(u* + e) = Tj. Most of the fixed-point 
subspaces in E are empty. The subgroups {Tj | Fix(Fj,.E) 7^ 0} can be arranged in a lattice of 
isotropy subgroups. The Equivariant Branching Lemma (EBL), described in [9], states that there 
is generically an EBL branch of bifurcating solutions with symmetry Tj if the fixed-point subspace 
Fix(Tj,E) is one-dimensional. In a gradient system such as PDE ([I]), there is generically a branch 
of bifurcating solutions with symmetry Tj if Tj is a maximal isotropy subgroup. See [9j [13] for 
details. 

The bifurcation digraph, defined in [13], summarizes some information about all of the generic 
bifurcations that are possible for a system with a given symmetry. In particular, if there is a 
daughter with symmetry Tj created at a generic bifurcation of a mother solution with symmetry 
Ti in a gradient system, then there is an arrow in the bifurcation digraph 

ryr' 

[r,]— 4[r,]. 

The arrow in the bifurcation digraph is either solid, dashed, or dotted, as described in [13]. Roughly 
speaking, a solid arrow indicates a pitchfork bifurcation within some one-dimensional subspace of 
E, a dashed arrow indicates a transcritical bifurcation within some one-dimensional subspace of E, 
and a dotted arrow indicates a more exotic bifurcation. 

In [13] we defined an anomalous invariant subspace (AIS) A C V, with V = Gm or H, to be a 
^-invariant subspace that is not a fixed-point subspace. Consider the PDE ([Tj) on the cube, with 
f s odd. For positive integers p, q, and r, not all 1, the set 

(11) A,q,r = span {^fc I |, j, J € Z} 

is an AIS. The function space -4i,i,i is all of V , so it is not an AIS. 

The space A P: q :T is the appropriate function space if one solved the PDE ([1]) in the box (0, ft/p) x 
(0, n/q) x (0, 7r/r). Any solution in this box extends to a solution on the cube, with nodal planes 
dividing the cube into pqr boxes. These anomalous invariant subspaces are caused by the so-called 
hidden symmetry [W\ [16] of the problem that is related to the symmetry of the PDE on all of M 3 . 

There are in fact a multitude of AIS for the PDE on the cube. There are proper subspaces of 
every AIS Ap^^ consisting of functions with symmetry on the domain (0, vr/p) x (0, ir/q) x (0, n/r). 
For example, consider the case where p = q = r > 1. Since there are 323 fixed-point subspaces for 
functions on the cube, including the zero-dimensional {0}, there are 322 AIS that are subspaces 
of A P)P)P , for each p > 1. Our code is capable of finding solutions in any AIS, within reason. The 
theory of anomalous solutions within AIS is unknown. The book [7] is a good reference on invariant 
spaces of nonlinear operators in general. 

2.3. Isotypic Decomposition. To analyze the bifurcations of a branch of solutions with symme- 
try Ti, we need to understand the isotypic decomposition of the action of Ti on R ra . 

Suppose a finite group T acts onV = R n according to the representation g i-> a g : T — > Aut(V) = 
GL„(K). In our applications we choose r 6 Q and the group action is the one in Equation (jlOj) . 

Let {a^ : T GL^^R) | k € K?} be the set of irreducible representations of T over R, 

where d r (k) is the dimension of the representation. We write and K when the subscript T 
is understood. It is a standard result of representation theory that there is an orthonormal basis 

= U k £K B r ] for v such that B r ] = Uz=i B ?' l) and [%U*.oL(W) = a(k) (g) for all g e T, 

where V^ k '^ := span(B^'^). Each V^ k ' 1 ^ is an irreducible subspace of V. Note that might be 

empty for some k, corresponding to = {0}. The isotypic decomposition of V under the action 
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of r is 

(12) ^=©^*\ 

where Vp = (Bz^i V^ k '^ are the isotypic components. 

The isotypic decomposition of V under the action of each Tj is required by our algorithm. The 
decomposition under the action of Aut(G) is the same as the decomposition under the action of 
To- While there are twice as many irreducible representations of Tq = Aut(G) x Z2 as there are 
of Aut(G), if apj(— 1) = / then VpT = {0}. The other half of the irreducible representations 
have a^\— 1) = —I. The irreducible representations of To and of Aut(G) can be labeled so that 

The isotypic components are uniquely determined, but the decomposition into irreducible spaces 

(k) (k) (k) 

is not. Our goal is to find l? r for all k by finding the projection P r : V — > Vp 1 . To do this, 
we first need to introduce representations over the complex numbers C for two reasons. First, 
irreducible representations over C are better understood than those over R. Second, our GAP 
program uses the field C since irreducible representations over M. are not readily obtainable by 
GAP. 

There is a natural action of T on W := C n given by the representation g 1— > (3 g : V — > Aut(W) such 
that j3 g and a g have the same matrix representation. The isotypic decomposition W = tB fcg ^- W r 
is defined as above using the set {/3^ k > : T —> GL^)(C) [ k € Kr} of irreducible representations of 
r over C. 

The characters of the irreducible representation (3^ are X (9) := Tr/3^ fc ^(<7). The projection 
: W -4 is known to be 



(13) Q? = %[^Zx {k) {g)Pg- 



(k) (k) 

We are going to get the P r 's in terms of Q-p s. A general theory for constructing these projections 
can be found in [13], but here it is enough that if = X^ k \ then = \y and dp = 

whereas if X (k) + W>, then P^ k) = (qP + Q (k) ^j |y and d (k) = 2d (k) , for all k G K T . 

2.4. Basis processing. In this subsection, we describe how the package of programs from [13] is 
modified to generate the basis needed to approximate solutions to PDE ([1]) on the square or cube. 
In principle a brute force method is possible, wherein the very large graph S7jv is used. However, 
the GAP and Mathematica portions of the package in [13] cannot process such large graphs with 
current computer systems. To remedy this, we wrote specialized Mathematica basis-generation 
programs for the square and cube. These programs use the known eigenfunctions, ipij or ipij^ki 
together with the GAP output from the graphs in Figure [1] to generate the data needed by the 
GNGA program. 

It can be shown that for polynomial f s and the known eigenfunctions ip m defined in terms of sine 
functions on the square and cube, the midpoint numerical integration can be made exact, up to 
the arithmetic precision used in the computation. In particular, consider the case when f s is cubic 
and the eigenfunctions are ipij (for d = 2) or ipij t k (for d = 3). Let M be the desired number of 
sine frequencies in each dimension, that is, i,j, k < M. Then our midpoint numerical integration 
is exact if 2M + 1 gridpoints in each direction are used, giving a total of N = (2M + l) d total 
gridpoints. For example, for analyzing PDE ([T]) on Q = (0, 7r) 2 we used M = 30 sine frequencies 
and, for exact integration, N = (2 • 30 + l) 2 = 3721 gridpoints. Similarly, for the cube Q = (0, 7r) 3 , 
we used M = 15 sine frequencies and N = (2 • 15 + l) 3 = 29, 791 gridpoints. 



PDE ON THE CUBE 



9 




R90 Rl20 Rl80 



Figure 2. Rqq, R120, and Riso are matrices that rotate the cube centered at the 
origin by 90°, 120° and 180°, respectively. The dashed lines are the axes of rotation. 
The 48 element automorphism group of the cube is identified with C\. The ma- 
trix group Oh is generated by these three rotation matrices and — 13, the inversion 
through the center of the cube. 

The only input for the basis generation program is M, the desired number of sine frequencies in 
each dimension. We define the bases 

B M = {ij> id I 1 < i,j < M and A*,,- = i 2 + j 2 < (M + l) 2 + 1} 

for the square, and 

Bm = {ip i;jjk I 1 < i,j,k < M and \ i>jjk = i 2 + j 2 + k 2 < (M + l) 2 + 2} 

for the cube. These bases include all eigenfunctions with eigenvalues less than A A 7 /+1 1 and A A/+1 1 x 

respectively. This process gives M = 719 for the square, not M 2 = 900, and M = 1848 for the 
cube, as opposed to M 3 = 3375. 

The basis generation code also produces an automorphism file, which is an 8 x (for the square) 
or 48 x N (for the cube) matrix describing how the elements of B4 or O/i permute the vertices in 
Oat. This file is used by the GNGA program but not by the square and cube basis generation 
code. Instead, the action of D4 or on the eigenfunctions is achieved very efficiently using the 
standard 2D representation of D4 acting on the plane and the standard 3D representation of 
acting on M 3 . This requires us to define tpij and V'ij.fc f° r negative integers i,j,k in this way: 
ip-ij(x,y) = ipi,j(Tr - x,y) and ip-ij >k (x,y,z) = ipi,j,k(^ ~ x, y, z). 

The GNGA program does not use the basis Bm- Rather, it uses a basis spanning the same 
M-dimensional space, obtained from projections of the eigenfunctions in Bm onto the isotypic 
components identified by the GAP program. The output of the GAP program for the small graphs 
shown in FigureCQis used to achieve these projections. For the cube, the eigenspaces of the Laplacian 
for the eigenvalue Aj k are 1-dimensional if i = j = k, 3-dimensional if i = j < k or i < j = k, and 
6-dimensional if i < j < k. Each of these eigenspaces are separated into their isotypic components 
automatically by the basis generation programs, using the characters found by GAP for the graphs 
in Figure [1] Then the projections in Equation f)13|) are written as 2 x 2 or 3 x 3 matrices, where 
fig is the standard 2D or 3D representation of B4 or O^, respectively. With these changes, the 
construction of the bases proceeds as described in |13j . 

2.5. The symmetry of a cube. In this section we describe the various matrix groups related to 
the symmetry of a cube. We start with the symmetry group of the cube (—1, l) 3 . The matrix form 
of this symmetry group is the 48-element := (-R90, R120, Rim-, ~H) < GL3(M) where 





"0 


-1 


0" 




"0 





1" 




"0 


1 


" 


-R90 — 


1 








J -R120 — 


1 








, -R180 — 


1 
















1 







1 













-1 



as shown in Figure [2j and I3 is the 3x3 identity matrix. 
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O h 2* §4 x Z 2 O ^ §4 T d ^ §4 T h ^ A 4 x Z 2 T ^ A 4 



Figure 3. Visualization of Oh and some of its subgroups. The group orbit of a 
single arrow is drawn, suggesting a symmetric vector field. These figures show the 
symmetry of the Lyapunov-Schmidt reduced gradient g on 3-dimensional critical 
eigenspaces in bifurcations with symmetry occurring in PDE (TjQ) on the cube. The 
dotted lines show the intersections of reflection planes with the cube. The dots 
indicate 1-dimensional fixed-point subspaces of the action, which intersect the cube 
at a vertex, edge, or face. The EBL states that solution branches bifurcate in the 
direction of these 1-dimensional fixed-point subspaces in E. 

There is a slight complication in describing the symmetry group of the cube O3 := (0, 7r) 3 that is 
the domain of our PDE, since is not centered at the origin. The action of Aut(f23) = on ^3 
is given by matrix multiplication about the center d = 7r/2(l, 1, 1) of the cube. That is, for 7 € Oh 
and x E Q3, the action is defined as 

7 • x = d + 7(x — d). 
The action of Oh x Z 2 on the vector space of functions u : O3 — > M is given by 
{{ 1 ,P)-u){x) = pu{ 1 - l -x), for 7 6 0/,, € Z 2 = {1, -1}. 
It follows that the action of the generators of Oh x Z 2 on eigenfunctions is 

(#90, 1) • 4>i,j,k = 
(#120 > 1) • ^i,j,k = i>h,i,j 
(#180, 1) • lfjij,k = (-l) fc ~Vj,i,fc 

(-/3,l)-V„,fc = (-ir +J+fc - 1 V'M,fc 
(I 3 , -1) • ifjij^k = -1pi,j,k- 

There are three subgroups of with 24 elements. They are 

O = (#90) #120 1 #180 ) ; 1^ = (— #90, #120) — #18o)) and Th = (#90, #120) #180) — ^3) 

shown in Figure [3l The group O contains the rotational symmetries of the cube (or octahedron) , 
and is called the octahedral group. Note that Oh is the internal direct product O x (—13). 
The groups and are related to the 12-element tetrahedral group 

= (#90> #120) #18o)) 

that contains all of the rotational symmetries of the tetrahedron. Note that is the internal direct 
product T x (—I3), whereas does not contain —I3. The groups O and are both isomorphic to 
the symmetric group S4. The isomorphism can be proved by considering how the groups permute 
the 4 diagonals of the cube. Similarly, T is isomorphic to A4, the alternating group. The group 
names involving §4 and A4 are used by the GAP program. In our GAP programs O and are 
computed as irreducible representations of §4. Our website, which is written automatically using 
the GAP output, also uses these names. 
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3. Algorithms 

The main mathematical algorithms used in the code for the results found in this paper are the 
tGNGA, cGNGA, and the secant method with recursive bisection, all of which are developed and 
described in detail in [13J. We give a brief overview of these algorithms in Section \3. 11 The current 
implementation of C++ code which supervises the execution of these algorithms has two substantial 
modifications not found in [13) . 

The first major improvement concerns the efficient way we now use our symmetry information 
to reduce the matrix dimension when setting up the system for the search direction x used in the 
tGNGA. In brief, each system uses only the rows and columns corresponding to the eigenfunctions 
having the symmetry of points on the branch. The time savings can be substantial when seeking 
solutions with a lot of symmetry, since the numerical integrations required to form the systems are 
generally the most time- intensive computations we make. We present the details in Subsection l3.2i 

Secondly, the current implementation is in parallel. A serial implementation could not reproduce 
our results in a reasonable time. To that end, in [15] we developed a simple and easy to apply 
methodology for using high-level, self-submitting parallel job queues in an MPI (Message Passing 
Interface [6]) environment. In that paper, we apply our parallel job queue techniques toward solving 
computational combinatorics problems, as well as provide the necessary details for implementing 
our PDE algorithms in parallel C++ code in order to obtain the results found in the current article. 
We include a high-level description in Section 13.31 

3.1. GNGA. To follow branches and find bifurcations we take the parameter s to be the (M + l) st 
unknown. When we say that p = (a, s) G R M+1 is a solution we mean that u = Yl m =i a m^ m solves 
Equation (j2]) with parameter s, that is g = 0. The (M + l) st equation, n(a,s) = 0, is chosen in 
two different ways, depending on whether we are implementing the tangent-augmented Newton's 
method (tGNGA) to force the following of a tangent of a bifurcation curve or we are implementing 
the cylinder-augmented Newton's method (cGNGA) to force the switching to a new branch at a 
bifurcation point. In either case, the iteration we use is: 

• compute the constraint k, gradient vector g := g s (u), and Hessian matrix h := h s (u) 



• solve 







Xa 
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X 3 . 




K 



• (a, s) <- (a, s) - x, u = Yl a j^r 
Equations © and ([8]) are used to compute g and h. The (M + l) st row of the matrix is defined 
by (VqKj = Vk G R M+1 ; the search direction is x = (Xa,Xs) £ K w+1 . Since this is Newton's 
method on (g, k) G M a/+1 instead of just g G 1 M , when the process converges we have not only 
that g = (hence p = (a, s) is an approximate solution to Equation ([T])), but also that n = 0. 

We use the tGNGA to follow branches. The details are given in Algorithm 1 of [13]. In brief, 
given consecutive old and current solutions p id and p CUT along a symmetry invariant branch, we 
compute the (approximate) tangent vector v = (p cur — p \d)/\\Pcuv — Poldll £ R + . The initial 
guess is then p gs = p CUT + cv. The speed c has a minimum and maximum range, for example from 
0.01 to 0.4, and is modified dynamically according to various heuristics. For example, this speed 
is decreased when the previous tGNGA call failed or the curvature of the branch is large, and 
is increased toward an allowed maximum otherwise. For the tGNGA, the constraint is that each 
iterate p = (a, s) must lie on the hyperplane passing through the initial guess p gs , perpendicular to 
v. That is, n(a, s) := (p — p gs ) ■ v. Easily, one sees that (V a n(a, s), s)) = v. In general, if f s 
has the form f s (u) = su + H(u), then || = —a. Our function tGNGA(p gs , v) returns, if successful, 
a new solution p new satisfying the constraint. Figure d] shows how repeated tGNGA calls are made 
when a worker executes a branch following job. 

The constraint used by the cGNGA (see Algorithm 3 of [13] ) at a bifurcation point p* instead forces 
the new solution p ncw to have a non-zero projection onto a subspace E of the critical eigenspace 
E. To ensure that we find the mother solution rather than a daughter, we insist that the Newton 
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repeat 

wait for a message from the boss 
switch the message is a 
case follow branch job 

while branch is in window do 

compute next point on branch with tGNGA 
if change in MI then 

call secant-bisection to find intervening bifurcation points 
for each intervening bifurcation point do 
|_ put find daughters job on queue 

]_ use interpolated guess and tGNGA to get last point on window boundary 

case find daughters job 

for each possible bifurcation subspace of critical eigenspace do 
while new solutions still being found do 

make a random guess in subspace and call cGNGA 
if new solution is nonconjugate to previously found solutions then 
[_ put follow branch job on queue 

case stop command 
|_ stop 



Figure 4. Pseudo code for the main loop of the workers. This loop is 
entered after loading basis and symmetry files. Whenever idle, each worker 
accepts and runs jobs whenever such jobs exist. The boss puts the trivial 
solution branch on the job queue as the first job. It manages the queue while 
the workers do their jobs, until the queue is empty, and then sends all workers 
a stop job. 

iterates belong to the cylinder C := {(a, s) G R M+1 : ||Pj;(a — a*)|| = e}, where Pe is the orthogonal 
projection onto E and the radius e is a small fixed parameter. At a symmetry breaking bifurcation 
the critical eigenspace is orthogonal to the fixed-point subspace of the mother, so the mother branch 
does not intersect the cylinder. The constraint we use to put each Newton iterate on the cylinder 
is K(a,s) = ^{\\Pe{o. — a*)\\ 2 — e 2 ) = 0. The initial guess we use is p gs := (a*,s*) +e(e,0), where 
e is a randomly chosen unit vector in E. Clearly, p gs lies on the cylinder C. A computation shows 
that V a K(a, s) = Pe(cl — a*), and §f («, s) = 0. When successful, cGNGA(p*, j> gs , E) returns a new 
solution Pncw that lies on the cylinder C. 

In the above paragraph, we take E to be various low-dimensional subspaces of the critical 
eigenspace, corresponding to the symmetries of solutions that are predicted by bifurcation theory. 
For example, at an EBL bifurcation E is spanned by a single eigenvector. When the dimension of E 
is greater than one, we call cGNGA repeatedly with several random choices of the critical eigenvector 
e. The details are given in Equation (7) and Algorithm 3 of [13]. The theory we apply does not 
guarantee a complete prediction of all daughter solutions. Therefore we also call cGNGA with E 
equal to the full critical eigenspace. In this way, if the dimension of the critical eigenspace is not 
too big we have a high degree of confidence that we are capturing all relevant solutions, including 
those that arise due to accidental degeneracy and that are neither predicted nor ruled out by un- 
derstood bifurcation theory. The number of guesses in each subspace is heuristically dependent on 
the dimension of the subspace. Too many guesses wastes time, and too few will cause bifurcating 
branches to be missed. Figure 0] shows how repeated cGNGA calls are made when a worker executes 
a find daughters job. 
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We use the secant method to find bifurcation points. In brief, when using the tGNGA to follow 
a solution branch and the MI changes at consecutively found solutions, say from k at the solution 
p \d to k + 5 at the solution p cur , we know by the continuity of D 2 J S that there exists a third, nearby 
solution p* where h is not invertible and the r th eigenvalue of h is zero, where r = k + [~| ] . Let 
Po = Po\d, Pi = Pcm, with j3o and /3i the r th eigenvalues of h at the points po and p\, respectively. 

We effectively employ the vector secant method by iterating 

_ (pj - pj-i)Pi 
•P&-Pi Uk _ fk _ l) 

• p i+1 = tGNGA(p gs , v) 

until the sequence (pi) converges. The vector v = (p C ur— Poid)/||Pcur— Poidll is held fixed throughout, 
while the value /3j is the newly computed r th eigenvalue of h at pi. If our function secant(p id,p C ur) 
is successful, it returns a solution point p* = (a*,s*), lying between p Q \d and p cnr , where h has 5 
zero eigenvalues within some tolerance. We take the critical eigenspace E to be the span of the 
corresponding eigenvectors. If p* is not a turning point, then it is a bifurcation point. 

In fact, it is possible that several intervening bifurcation points exist. If the secant method 
finds a bifurcation point that has fewer than d zero Hessian eigenvalues, there must be another 
bifurcation point in the interval. In Figured! the secant-bisection call refers to an implementation 
of Algorithm 2 from p.3] entitled f incLbif points, whereby such an occurrence triggers a bisection 
and a pair of recursive calls back to itself. Upon returning, each of the one or more found bifurcation 
points spawns its own find daughters job. In turn, each time a (non-conjugate) daughter is found, 
a new follow branch job is put on the queue. 

3.2. The Block Diagonal Structure of the Hessian. The majority of the computational effort 
for solving PDE ([1]) using Newton's method comes from the computation of the entries of the 
Hessian matrix. The time required can be drastically reduced by taking advantage of the block 
diagonal structure of the Hessian that follows from the isotypic decomposition of V = 1R M , the 
Galerkin space. 

If the initial guess u has symmetry Tj, then the isotypic decomposition ()12|) of the Tj action on 
V is V = ® fc Vp. , where k labels the irreducible representations of Tj. We assume that k = 

denotes the trivial representation, so Fix(Fj) = Vp . 

For any u G Fix(Fj), the symmetry of the PDE implies that the gradient g s (u) is also in Fix(Fj), 
and the Hessian, evaluated at u, maps each of the isotypic components to itself. That is, 

u G Fix(Ti) g s (u) G Fix(Ti) and h s (u) (v£A C vff. 

Thus, the Hessian is block diagonal in the basis Br l defined in Section 12.31 A huge speedup of 
our program is obtained by only computing the Hessian restricted to Fix(Tj) when doing Newton's 
method. After a solution is found, the block diagonal structure of the Hessian allows its efficient 
computation by avoiding integration of zero terms. The full Hessian is required for the calculation 
of the MI. 

Actually, the way we achieve a speedup in our numerical algorithm is not quite this simple. In 
our implementation of the GNGA we always use -E>r > a basis of eigenvectors of the Laplacian. The 

(k) 

basis vectors are partitioned into bases Bf, for each of the isotypic components of the Tq action on 
V. We do not change the basis depending on the symmetry of the solution we are approximating. 
Hence, we do not simply compute the blocks of the block diagonal Hessian. 

When doing Newton's method, we use a reduced Hessian h s in place of the full Hessian h s . Define 
p( k ">v to be the projection of v onto vS k \ For u G Fix(Fj), the reduced Hessian is defined to be 



'h s {u)j, k if p(°Vj + and P<-°ty k + 
Xj -s if j = k and P^Vj = 
if j ^ k and {P&tyj = or P(°tyk = 0). 
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For each u € Fix(Tj), assuming h s and h s are nonsingular, the Newton search direction is the 
solution to either system h s (u)x = g s (u) or h s (u)x = g s (u). The terms in the reduced Hessian of 
the form (Aj — s)5jk are included to make h s nonsingular. They are not strictly necessary since we 
find the least square solution \ with the smallest norm, but they improve the performance of the 
LAPACK solver (dgelss). 

After the solution is found, we identify those elements of the full Hessian that are known to be 
zero, and avoid doing numerical integration for those elements. In particular, h s (u) p ^ q = if there 
is no k such that ^ and P^^j q ^ 0. 

As a simple example for demonstration purposes, suppose that there are M = 5 modes, and we 
are using the basis vectors ipi, . . . ,^5. Suppose further that the isotypic decomposition of the r. 



action is Vp = span{^i, ^3}, Vp X} = span{-02}, V^' = span{^4 + V's}) an d Vp°' = span{^4 — ^5}. 
The only nonzero projections are P^ipi, P^ip s , P^ip 2 , P^ip^, P^ips, P^ip^, and P^ip 5 . Thus, 
if u G Fix(Tj), the gradient and Hessian have the form 



r(l) 



(2) 



,(3) 



(14) 



9s{u) 



h s (u) 



* * 
0*000 

* * 
** 
** 
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A 2 - s 
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* 

















A4 — s 

















A5 - s 



Note that a change of basis could be done to diagonalize the lower right 2x2 block, but our 
program does not do this. When we solve for the search direction % in Newton's method, we use 
the restricted Hessian 



h s (u) 



where only the * terms are computed using numerical integration. Thus, in this small example, 
four numerical integrations are needed to compute the reduced Hessian for each step of Newton's 
method, and nine numerical integrations are needed to compute the full Hessian. 

The speedup obtained in this manner is quite dramatic for solutions with high symmetry when 
M is large. As an example, for the PDE on the cube the primary branch that bifurcates at s = 3 has 
symmetry T2 — Oh- The T2 action on V has 10 isotypic components. Our program, with M = 15 
and hence M = 1848, processes these solutions 10 to 15 times faster than it processes solutions 
with trivial symmetry, where there is just one isotypic component and all M(M + 1)/2 = 1, 708, 474 
upper-triangular elements of the Hessian need to be computed for every step of Newton's method. 
Each one of these Hessian elements requires a sum over = 29, 729 grid points (see Equation ([!])). 
Considering the whole process of finding a solution with T2 symmetry using 3 or 4 steps of Newton's 
method, and then computing the MI of this solution, the majority of the computational effort comes 
from the numerical integrations needed to compute the full Hessian once after Newton's method 
converges. Even with the speedup indicated in (|14h . the single computation of h s (u) takes longer 
than the total time to perform all the other computations. These include the computation of the 
reduced Hessian h s (u) and gradient g s (u), and the LAPACK calls to solve h s (u)x = g s (u) at each 
Newton step, as well as the LAPACK computation of the eigenvalues of the full Hessian evaluated 
at the solution. 



3.3. Parallel Implementation of Branch Following. It was necessary to implement our code 
in parallel in order to get timely results. We use a library of functions (MP Queue) presented in [15] 
which make it easy to create and manage a parallel job queue. For the current application of 
creating a bifurcation diagram, we choose a natural way to decompose the task into two types of 
jobs, namely branch following and find daughters jobs. The boss starts MPI and puts the trivial 
solution branch on the queue as the first job. After initialization and whenever idle, each worker 
accepts and runs jobs whenever they exist. The boss manages the queue while the workers do their 
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Figure 5. The load diagram showing the number of jobs (branch following and 
bifurcation processing) and the number of active workers. There were 24 processors 
and thus 23 workers. The jobs above the dotted line are waiting in the job queue. 
The total time for the run was 192,094 seconds ~ 53.4 hours. The average number 
of active workers is 15.6. Thus, the run would have taken approximately 34 days 
with a single processor. 



jobs, until the queue is empty, and then sends all workers a stop command. Figure U] shows how a 
parallel job queue is used to supervise the running of the jobs. 

Normal termination for branch following jobs occurs when the branch exits a parameter interval 
for s. When a change in MI is detected during branch following, the secant method is used to find 
an intervening bifurcation point with a proscribed zero eigenvalue of the Hessian of J. Recursive 
bisection is used if the number of zero eigenvalues of a found bifurcation point does not equal the 
observed change in MI over a given subinterval, ensuring that all intervening bifurcation points 
are found. Each bifurcation point spawns a find daughters job. These jobs invoke the cGNGA 
with multiple random guesses from all possible bifurcation subspaces of the critical eigenspace. 
The number of attempts depends heuristically on the dimension of each subspace, stopping after 
sufficiently many tries have failed to find a new solution not in the group orbit of any previously 
found solution. Each new, bifurcating solution spawns a follow branch job. 

We conclude this subsection with an example demonstrating the parallel run times encountered 
when generating our results for the cube region. In particular, Figure [5] shows the load diagram 
of a run for PDE ([!]) on the cube with M = 15, which gives M = 1848 modes and uses N = 31 3 
gridpoints. We followed all branches connected to the trivial branch with < s < 13, up to 4 
levels of bifurcation. The code found 3168 solutions lying on 111 non-conjugate branches, with 126 
bifurcation (or fold) points. This 2-day run used 24 processors and generated all the data necessary 
for Figures [1411151 1201 1211 and [22] In Table 6.1 of [15] we provide runtime data for a similar problem 
using a varying number of processors, in order to demonstrate scalability. 

3.4. Creating contour plots. We now describe the process of creating a contour plot for visually 
depicting an approximate solution on the square and cube. For the square contour plot graphics 
found in Figures [9] and [Til we use M = 30 leading to N = 61 2 grid points. This is a sufficient 
number of grid points to generate good contour plots. The cube solutions visualized in the contour 
plot graphics found in Figures [151 ttH GU [2TJ and [25] were all obtained via GNGA using M = 15 
leading to M = 1848 modes and N = 31 3 grid points. This resolution is not quite sufficient for 
generating clear contour plots. Using the M coefficients and the exact known basis functions, we 
reconstruct u on a grid with 56 3 points before applying the techniques described below. 
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( e ) -w (f) -w (g) 

Figure 7. Boundary detection, followed by polygon creation and polygon merging. 



When Q is the cube, since u = on the boundary we plot u on the boundary of a smaller 
concentric cube. The faces of this inner cube are 5 gridpoints deep from the original boundary. 
Thus each of the 6 faces we draw has u values on (56 — 2 • 5) 2 grid points. The side length of the 
cube seen in the figures is about 46/56 = 82% the length of the original cube. At the end of the 
process, the contour plots on the six faces of the inner cube are assembled into a front and back 
view. We use a standard perspective mapping from M 3 to R 2 to join three of the six faces for each 
view. 

In the remainder of this sub-section, we describe how to make a contour map for a function 
u : ft — > K on the square or cube. First of all, an algorithm chooses one function in the group 
orbit Tq ■ u that makes the symmetry of the contour map agree with the flag diagram. It also 
replaces u by — u if this saves ink. For example, if the function has a 3-fold symmetry about a 
diagonal of the cube, our algorithm chooses a function that satisfies (Ri20> l)-u = u since the axis of 
rotation of R120 points toward the viewer. Then we have the values of u on a square grid (from the 
square or a face of the cube) as shown in Figure [6(a) . We do some data pre-processing, depicted in 
Figures [6(b-d) , in order to get acceptable plots when there are many data points with u = 0, which 
is common because of the symmetry of solutions. After interpolation, Figure [7(e), the positive 
region is the union of many polygons, mostly squares, Figure [7(f). Redundant edges are eliminated 
to obtain a few many-sided polygons that are shaded, Figure [7(g). The contours are obtained 
by the same algorithm: It finds the many-sided polygons enclosing the regions where u > c, but 
just the boundary is drawn. Our implementation directly creates postscript files. Without the 
elimination of redundant edges, the files would be of an unreasonably large size, and the contours 
could not be drawn with the same algorithm. 
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We now provide the details. Let (x, a) represent a point where x £ M. n , n = 2 or 3, is a grid 
point and a£l such that u(x) ~ a for a given solution u. We use the following steps as visualized 
by Figures [6] and [71 

(1) Square creation (a)^»(b): In the first step, we break the square grid into individual 
squares to be handled separately. We represent one of these squares with the cycle 
((xi, ai) • • • (X4, 04)) containing the corner points. 

(2) Zero localization (b)-w(c): We eliminate several zero values in a row by applying the rule 

• • • (x,a)(yi,0) • • • (y n ,0)(z,6) (x, a) (yi , a) • • • (y„,6)(z,6) • • • 

to each cycle where the variables a and b represent non-zero values. 

(3) Zero purge (c)~»(d): We eliminate the zero values without a sign change by applying the 
rule 

• • • (x, a)(y, 0)(z, b) >■■■ (x, a)(y, (a + 6)/2)(z, b) ■ ■ ■ 

where ab > 0. 

(4) Boundary detection (d)-w(e): We insert zero points between sign changes using the rule 

. . . ( X , a )(y,6) ....... (x,o) ( |6 |^+ j°j y , 0) (y,6) ■ • • 

where ab < 0. The location of the new zero point is determined using linear interpolation. 
The number k of zeros in the cycle after this step must be 0, 2 or 4. 

(5) Polygon creation (e)~»(f): Now we create polygons with shaded inside to indicate positive 
values. 

(a) If k = and ai > for all i in the cycle ((xi,ai) • • • (x4,04)), then we create the 
polygon with these corner points. This occurs in the lower-right square of Figure 0(e) . 
Note that if a, < then no polygon is created. 

(b) If k = 2, then we create a polygon for each cyclic pattern of the form 

(x,0)(yi,ai) • • • (y n ,a n )(z,0) 

with ax, . . . , a n all positive. There is an example of this with n = 2 in the upper-right 
square and an example with n = 3 in the lower- left square of Figure [TJe) . It is also 
possible to have n = 1, but this is not shown in the example. 

(c) If k = 4, then we find the intersection (c, 0) of the line segments joining opposite zeros. 
Then we look for cyclic patterns of the form 

(x,0)(y,a)(z,0) 

with a positive. For each such pattern we create the polygon with corner points 

(x,0),(y,a),(z,0),(c,0). 

This happens in the upper-left square of Figure [T^e) 

(6) Polygon merging (f)~»(g): We combine polygons sharing a side into a single polygon. This 
gives the zero set which is the boundary of the shaded positive region. This single polygon 
is written directly to a postscript file. Polygon merging is a time consuming operation, but 
writing the merged polygons instead of the individual ones reduces the size of the postscript 
file significantly. 

(7) Level curve creation: Level curves are created using the previous steps. We use the zero set 
of u — c as the level curve of u = c. The original grid values are changed by the shift value 
c and only the final merged polygon is written to the postscript file without shading. Note 
that polygon merging is essential for finding level curves. 

(8) Local extremum dots (g): White dots are drawn at the estimated position of local extrema 
in the u > region, as indicated in Figure E^g). Every interior data point (x, a), is checked 
to see if a is at least as large (or at least as small) as 6, for the eight neighbors (xj,6j). For 
the data in Figure 0(a), only the center point is checked, and it is indeed a local maximum 
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grid point. For each local extremum grid point, a more precise estimate of the extreme 
point is found by fitting a quadratic function to the 9 data points with the extreme grid 
point at the center. Since a quadratic function / : M 2 — > M has 6 constants, the least 
squares solution for the 9 equations in 6 unknowns is found. Then standard calculus gives 
the position of the extreme point of /, and a dot is drawn there. For the data in[6{a), 
the maximum point of the fitted / is located half way to the grid point to its lower left. 
Similarly, black dots are drawn at the estimated positions of the local extrema in the white 
(u < 0) region. 

4. The PDE on the Square and B 4 Symmetry 

In this section we re- visit the case where O is the unit square. We first studied the square case 
in |16j . without the benefit of automation or our recent improvements in branch following. In that 
paper, all of the symmetry analysis was done by hand. Each of the few bifurcation points we 
analyzed were essentially small projects in themselves, as were monotonic segments of the branches 
connecting them. In the present article, we use our automatically generated bifurcation digraph and 
new Newton algorithms to quickly reproduce the old results and then go much further (s S> Ai) in 
finding solutions of every possible symmetry type. We also show that we can handle some accidental 
degeneracies that occur in the square case. 

Figure [8] contains the condensed bifurcation digraph for PDE ([1]) on a region with D4 symmetry. 
This digraph, along with the files required by our GNGA code for following bifurcations, was 
automatically generated by our code. For clarity, we have chosen here to annotate the vertices 
of the digraph with schematic diagrams, i.e., contour plots of step functions on the square which 
display the proper symmetries in a visually obvious way. A solid line in the schematic diagrams 
represents a nodal line, whereas a dashed line is a line of reflectional symmetry and a dot is a 
center of rotational symmetry. Contour plots of actual solutions can be found in Figure [H The 
bifurcation digraph was described in Subsection 12.21 and a more thorough discussion is in |13| . 

The bifurcation digraph in Figure [8] is condensed, as described in |13| . The symmetry types 
are grouped into condensation classes, and not all of the arrows are drawn. For example, one 
condensation class is the block of four symmetry types near the top. The condensation class has 
three arrows emanating from it, but in the un-condensed bifurcation digraph there are 5 arrows 
emanating from each of the 4 solution types in the block. The little numbers near the arrow tails 
count the number of arrows emanating from each solution. Similarly, the little numbers near the 
arrow heads count the number of arrows ending at each of the solution types in the condensation 
class. 

At the top of Figure [8] is the trivial function u = whose symmetry is all of B4 x Z2. At the 
bottom is a function with trivial symmetry, i.e., whose symmetry is the group containing only 
the identity. There are four generic bifurcations with Z2 symmetry from the trivial branch. For 
these bifurcations the critical eigenspace E is the one-dimensional irreducible subspace for Z2 and 
there is a pitchfork bifurcation creating two solution branches, (u s ,s) and (— u s ,s) at some point 
(0,s*). The figure at the bottom left indicates the symmetry of a vector field in E. We can 
think of this as the ODE on the one-dimensional center manifold, or Lyapunov-Schmidt reduced 
bifurcation equation g = 0. There is one generic bifurcation with D4 symmetry from the trivial 
branch. Here the two-dimensional critical eigenspace E has lines of reflection symmetry across 
"edges" and non-conjugate lines of reflection symmetry across "vertices" , as shown in the bottom 
middle part of Figure At the bifurcation there are two conjugacy classes of solution branches, 
therefore the bifurcation digraph has two arrows labeled D4 coming out of the trivial solution. On 
the condensed bifurcation digraph these two arrows are collapsed into one. Note that there are 
several more bifurcations with Z2 symmetry or with B4 symmetry in the bifurcation digraph. For 
example, each of the 4 solution types in the second row can have a bifurcation with D4 symmetry. 

There is only one more type of generic bifurcation that occurs in PDE ([1]) on the square: a 
bifurcation with Z4 symmetry. This is a generic bifurcation in this gradient system |13j . and 




Figure 8. Condensed bifurcation digraph for PDE on the square, and the irre- 
ducible spaces for the generic bifurcations in the digraph. 




Figure 9. Contour plots of solutions on the square of each symmetry type, with 
schematic diagrams (see Figure [8|). 
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Figure 10. A partial bifurcation diagram for PDE ([I]) on the square, showing 
just one primary branch and selected daughter branches. The value of u(x*,y*) vs. 
s is plotted, where (x* , y*) is a generic point of the square, as shown in the figure on 
the right. A generic point is not on any of the lines of reflection symmetry, and we 
choose a point equidistant from those lines and the boundary. The open dots show 
bifurcation points, and the group of the bifurcation is shown for some. The small 
numbers indicate the Morse Index of the solutions, and solid lines show an even MI 
whereas dashed lines show an odd MI. The solid dots at s = 0, from top to bottom, 
correspond to the contour plots in Figure HH from left to right. 

the daughter solutions can be anywhere in the two-dimensional irreducible subspace (except the 
origin) since there are no lines of reflection symmetry. Note the use of a dotted arrow type for this 
bifurcation, and that there are no dashed arrows in this particular bifurcation digraph. The reduced 
bifurcation equations on the critical eigenspace at this bifurcation has the symmetry indicated at 
the lower right part of Figure 

Figure [9] contains the contour plot of an example solution to Equation (pQ) at s = for each of 
the 20 possible symmetry types on the square. The contour heights are ±c2~ h with h £ {0, . . . , 4} 
and an appropriate c near max(|ii|), to give more contours near u = 0. These figures were made 
with M = 30, meaning that the largest frequency in each direction is 30. Thus the mode with the 
smallest eigenvalue that is left out of the basis is ^31,1- This leaves M = 719 modes in our basis. 
Unlike [16], where an initial guess for each branch needed to be input by humans, the solutions in 
Figure [9] were found automatically by following all of the primary branches that bifurcated from 
the trivial solution, and all of the secondary branches, etc., recursively as described in |13j . 

In addition to generic bifurcations, the PDE on the square has degenerate bifurcations due to 
the "hidden symmetry" of translation in the space of periodic functions [10]. For example, the 
bifurcation point at s = X^^ = 34 on the trivial branch a* = has an accidental degeneracy of 
Type 2 as defined in [13]. Figure [TOl shows a partial bifurcation diagram containing this point 
and 3 levels of branches bifurcating from it; corresponding contour plots of solutions are found in 
Figure [TTJ A bifurcation diagram showing branches that bifurcate at s < A2 3 = 13 is shown in 

The 2-dimensional critical eigenspace at this primary bifurcation point is E = spanj^s, ^5,3}- 
The trivial branch, whose symmetry is Tq = B4 x Z2, undergoes a bifurcation with Tq/T' q = Z2 x Z2 
symmetry at s = 34. The action of Z2 x Z2 on E is generated by 



¥3,5 + ^5,3^^5 + 6^5,3 and 6^3,5 + c^5,3 ^ -^3,5 - ¥5,3 • 
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Figure 11. A sequence of solutions obtained at s = by following a chain of 
Z2 bifurcations from the trivial solution with full D4 x Z2 symmetry, down to a 
solution with trivial symmetry. The sequence of bifurcations shown is a path in the 
bifurcation digraph in Figure O The first solution shown is on a primary branch 
bifurcating at s = Xs^ = = 34, and is approximately a multiple of ^3,5 — ^5,3 • 
The four contour plots represent the solutions indicated by the four black dots in 
Figure [TOJ 

Our code uses the ordered basis (^3,5 + ^5,3,^3,5 — ^5,3)- in this basis, the action of Tq/F' on E 
is isomorphic to the natural action of 

([o- i] 5 [V-°i]) = ^xZ 2 

on [E] = IR 2 = K © R. Note that E is not an irreducible space; this is a degenerate bifurcation. 
The basis vectors were chosen to span the two one-dimensional irreducible subspaces, which are 
also fixed-point subspaces of the Tq/T'q action on E and therefore ^-invariant subspaces. Each of 
these subspaces is used as an E C E in the cGNGA algorithm described in Section 13. 1L A pitchfork 
bifurcation generically occurs in each of these one-dimensional invariant spaces. 

The one-dimensional subspaces spanj^s} and span{^5,3} Q E, corresponding to span{(l,l)} 
and span{(— 1,1)} C M 2 , respectively, are not fixed-point subspaces. However, they are AIS (see 
Section I2.2[) of g : E — > E, since V>3,5 (and ^5,3) can be periodically extended to tile the plane 
with a solution to the PDE ([1]). There is a primary branch of solutions which is tangent to 
{(u,s) = (0^3,5, 34) a £ t} at (0,34). Thus, there are at least three (conjugacy classes of) 
solution branches bifurcating from this degenerate bifurcation with Z2 X Z2 symmetry. Near the 
bifurcation, the nontrivial solutions are approximately multiples of ^3,5 + ^5,3 > ^3,5 — ^5,3 > or V>3,5 
(or its conjugate ^53). Figure [TOl shows a partial bifurcation diagram which follows the primary 
branch which is approximately a multiple of ^5 — ^5,3 near the bifurcation. 

Figure [TT] shows contour plots of example solutions along a particular path in the bifurcation 
digraph shown in Figure El The primary branch is created at the degenerate bifurcation with 
Z2 x Z2 symmetry at s = 34 discussed above. The critical eigenspace is two-dimensional at s = 
A3 t 5 = 34, and this bifurcation is not on the bifurcation digraph, Figure which only shows generic 
bifurcations. There are two primary branches conjugate to the one shown, four secondary branches, 
eight tertiary, and 16 branches conjugate to the solution with trivial symmetry shown in Figure HU 
At each bifurcation our GNGA code follows exactly one of the conjugate branches that bifurcate. 

Figure [12] contains a numerical demonstration of the convergence of the GNGA as the number 
of modes increases. In Figure 7(b) of [16j . we previously provided a portion of a similar graph 
of the L 2 norm of Au + u 3 vs. M, for a particular solution u. That graphic was not entirely 



PDE ON THE CUBE 



23 
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Figure 12. A measure of the error of our approximation to the CCN solution 
(minimal energy sign-changing solution [1]) of PDE ([1]) at s = 0, as a function of 
M. The region used was ft = (0, l) 2 to facilitate comparison with Figure 7(b) of 

M- 

convincing in showing convergence to 0. Using our current code with larger values of M and a 
smaller convergence tolerance in Newton's Method, we re-computed the data for the same problem 
and obtained the more accurate numerical result displayed here. 

5. The PDE on the Cube 

The PDE (fTJ on the cube has a rich array of bifurcations with symmetry. If fi is a planar region, 
only bifurcations with D n or Z n symmetry are present in the bifurcation digraph. The bifurcation 
digraph of the PDE on the cube includes bifurcations with Z2 symmetry, which have 1-dimensional 
critical eigenspaces. The digraph also includes bifurcations with Z n or D n symmetry, n € {3,4,6}, 
for which the critical eigenspace is 2-dimensional. A novel feature of the PDE on the cube is 
that there are 5 bifurcations with symmetry that have 3-dimensional critical eigenspaces, shown in 
Figure El 

Recall |13j that a generic bifurcation with T symmetry has a critical eigenspace E that is a faithful, 
irreducible representation space of T. Faithful means that only the identity in T acts trivially on 
E, and irreducible means that no proper subspace or E is T-invariant. 

Faithful, irreducible representation spaces for ID 4 and Z4 were shown in Figure El and the gener- 
alization to D n and Z n is obvious. Note that there are no such representation spaces for Z2 x Z2, 
since no 1-dimensional representation space is faithful, and every 2 or higher-dimensional represen- 
tation space is reducible. Recall that the bifurcation with Z2 x Z2 symmetry that occurs at s = 34 
in Figure [10] is not generic. That bifurcation point has a Type-2 degeneracy |13] . 

This section contains our main new numerical results, namely approximate solutions to Equa- 
tion (fTJ on the cube. For convenience, we denote the 99 symmetry types of solutions to this PDE on 
this region by fio, . . . , Sgs- Obtaining accurate approximations on this 3-dimensional region requires 
a large number of gridpoints. We use a parallel implementation as described in Section 13.31 

In Subsection 15. R we give an overview of features of the bifurcation digraph, which is too big to 
include in its entirety in a single document. We describe how our companion website [14] can be used 
to navigate the digraph in order to view graphics and understand various symmetry information 
across the spectrum of solutions. In the remaining subsections, we include a survey of our numerical 
results which showcase our analysis of the bifurcations with most interesting symmetries. 

5.1. The Bifurcation Digraph. The bifurcation digraph of the O/j x Z2 action on V = Gm or 
H is far too complicated to display as a figure in this paper. Our web site [13] has a page for each 
symmetry type Si, encoding the arrows emanating from this symmetry type, along with additional 
information. 
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Symmetry types, partitioned into condensation classes. 
Non- isomorphic types are separated with a ':' 

:0: 1 2 3 4 : ■ ■ ■ : 19 : 20 21 22 23 : ■ ■ ■ : 87 : 88 . 89 90 91 92 . 93 94 . 95 . 96 97 : 98 
Bifurcates from: 2 

Symmetry type: 21, representatives are isomorphic to Bg 
(Flag diagram and back and front contour plots appear here.) 
Bifurcation with (symmetry) to (symmetry type): 

| Z 2 ->• 55 | Z 2 ->■ 58 | Z 2 -» 57 | B 6 ->• 91 92 | B 3 > 80 

View: 55 58 57 91 92 80 



So 5 2 




Figure 13. A schematic representation of the page for symmetry type S21 from 
the companion web site |14j . along with the corresponding arrows in the bifurcation 
digraph. The labels on the top two arrows are found in the pages for symmetry 
types So and 5 2 , respectively. 

An example from the companion web site for symmetry type 5 2 i is shown in the top half of 
Figure [T3l The first box contains links to all the symmetry type pages. The 99 symmetry types are 
grouped into isomorphism classes by colons. The isomorphism classes are further subdivided into 
condensation classes by periods. The abbreviated list in Figure [TBI shows, for example, that So, S±g, 
Ss7 and 599 are singleton condensation classes, and that {620, S21, S22, $23} is a condensation class. 
The symmetry types S^s through S97 are all isomorphic. These 11 symmetry types are separated 
into 5 condensation classes. For example {5*93,594} is a condensation class. 

The second box indicates that there are arrows in the bifurcation digraph pointing from 5o and 
from S2 to S21. 

The third box indicates that any Tj € S21 is isomorphic to Oe> an d contains the graphics on the 
web page. There is a flag diagram for every symmetry type, and contour plots if our computer 
program found a solution with this symmetry type. About half of the symmetry types feature 
contour plots. 

The fourth box encodes the 6 arrows in the bifurcation digraph emanating from 5 2 i, as shown in 
the bottom half of Figure [T3j In addition, the arrows are separated into the 5 generic bifurcations 
with symmetry coming from the 5 nontrivial irreducible representations of B@. The symmetry types 
in any condensation class have an identical pattern of generic bifurcations, with different labels of 
the symmetry types. Thus, there are 6 arrows emanating from each of the symmetry types in the 
class {5 2 o, 5 2 i, 5 22 , 5 2 3}. In the condensed bifurcation digraph, these 24 arrows are represented by 
just 6 arrows. 

The final box contains buttons to view selected daughter flag diagrams. 

5.2. Bifurcations from the first three eigenvalues. Figure [141 shows the bifurcation diagram 
for the solution branches that are connected to the trivial solution branch with s < 10. The value 
of u at a generic point is plotted against the parameter s. The trivial branch has a bifurcation 
with Z 2 symmetry at s = \x t i,i = 3 where the MI changes from to 1. The trivial branch has two 
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Figure 14. Partial bifurcation diagram for PDE ([JJ on the cube, showing the 
first three bifurcation points of the trivial branch. The two bifurcations with Oh 
symmetry and three bifurcations with 7Li symmetry are indicated by open circles. 
The symmetry type of the trivial branch is Sq, as indicated on the right. The 
symmetry types of the other branches are shown on the left. The small numbers 
indicate the Morse Index of the solution. The solution branches with symmetry 
type 52i and S45 are truncated to simplify the diagram. The solid dots at the left 
correspond to the contour maps shown in Figure \15\ and contour maps for the other 
solution branches are shown in Figure [TBI 

bifurcations with Oh symmetry at s = Ai 1 2 = 6 and at s = Xi,2,2 = 9. The MI changes by 3 at 
each of these bifurcations, indicating a 3-dimensional critical eigenspace E. 

At s = 3, the symmetry of the mother branch is Tq = Oh X Z2. The critical eigenspace is 
E = spanj-^i^i}. All of the reflections and rotations in act trivially on E, That is, T' = 
((R90, 1), (-Ri20> 1)) (-^180) I); (-^3> 1)} = F2 — Oh- Thus, the effective symmetry of the bifurcation is 
Tq/Tq = ((I3, — l)r' ) = Z2 = { — 1,1}. The primary branch created at s = 3 has symmetry T2, and 
symmetry type S2 = ^2}, as shown in Figures [TH and [T5l 

The trivial solution undergoes a bifurcation with Oh symmetry at both s = \i,i2 = 6 and 
s = Xi t 2,2 = 9. The bifurcations are very similar. The critical eigenspaces for the two bifurcations 
are 

E 6 = span{V>2,i,i, ^1,2,1, ^1,1,2}, and E 9 = span{V>i,2,2, ^2,1,2, ^2,2,1}, 

respectively. The kernel of the action of Tq = Oh x Z2 on E% is T' = ((I3, —1)) and the kernel of 
the action of Fq on Eg is T' = ((—I3, —1)}- In both cases, 

ro/r'o = ((Rqo, i)r' , (r 120 , i)r' 0) (r 1B0 , i)V ) = o h . 

The reduced gradient g : E — > E in both cases has the equivariance indicated in the first image of 
Figure El The Equivariant Branching Lemma (EBL) guarantees that under certain non-degeneracy 
conditions there is a bifurcating branch tangent to each one-dimensional fixed-point spaces in E. 
These fixed-point subspaces intersect a cube in E at the center of a face, the center of an edge, 
or a vertex of the cube. Thus, each EBL branch is made up of face, edge, or vertex solutions. 
The standard choice of representative in each symmetry type is Tj € Si, for < i < 98. The 
contour plots of the bifurcating solutions in Figure PT51 show the solution in the fixed-point subspace 
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Figure 15. Contour plots for the solutions indicated with a solid dot in Figure HH 
The contour lines for functions on the cube are equally spaced. The positive solution, 
with symmetry type S2, bifurcates at s = 3. Representatives of the face, vertex, 
and edge solutions (denoted by F, V and E) that bifurcate at s = 6 are shown. The 
vertex solution with symmetry type S22 is the CCN solution. The names derive for 
the position of the critical points of g : E — > E. The equivariance of the vector field 
g is shown in Figure [3l 

indicated here: 

face [Fix(r 12 ,£ 6 )] = [Fix(r 14 , Eg)] = {(0,0, a) | a € R}, 

vertex [Fix(r 22 , E 6 )} = [Fix(r 2 i, Eg)] = {(a, a, a) | a G R}, 

edge [Fix(r44, Eq)] = {(— a, a, 0) | a G R} is conjugate to 

[Fix(r 45 ,^9)] = {(a, a,0) I a G R}. 
The symmetry type containing Ti 2 , denoted S*i 2 , has 3 elements, and the three conjugate face 
directions in [E] are the coordinate axes. Similarly, [F 22 ] = 5 22 has 4 elements, corresponding to 
the four diagonals through vertices of the cube centered in R 3 . The edge solutions bifurcating at 
s = 6 have symmetry type ^44] = £44, which has 6 elements. 
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FIGURE 16. Contour plots for the vertex, edge, and face solutions that bifurcate at 
s = 9 in Figure [HI along with a daughter and granddaughter of the face solution. 
The contour plots are for the solutions at s = 0. 
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Figure 17. Bifurcation with J}q symmetry in PDE JIJ on the cube, at s = 11. 
The symmetry of g on the critical eigenspace E is shown in Figure [181 There are 5 
(conjugacy classes of) branches that bifurcate at s = 11. One branch in each class 
is shown in this bifurcation diagram, and the solid dots indicate the solutions with 
the contour plots in Figure [T9l (The upper solid dot is very close to a bifurcation 
point.) There is an anomaly-breaking bifurcation at s ~ 3.417, where mother and 
daughter both have symmetry type S17. 



For the bifurcation at s = 6, the geometry in Eq is mirrored in the geometry of the solutions. 
For example, the maximum u value for a vertex solution (type S22) lies on the line from the origin 
to a vertex in 0,. Similarly, the maximum u value for an edge or face solution is on the line from the 
origin to an edge or face, respectively. The three face directions in E are ^2,1,1 j "01,2,1) an d ^1,1,2 > 
which can be thought of as u x" , "y", and u z" functions. Note that the face solution in Figure [15] 
is approximately a multiple of ^1,1,2 • 

The geometry in Eg for the bifurcation at s = 9 is the same. The geometry of the bifurcating 
solutions, shown in Figure [TBI is more subtle though. The face solutions have symmetry type S14, 
and the "z" eigenfunction is ^2,2,1 • The face and edge solutions have a line where two nodal planes 
intersect at right angles, and these lines intersect the midpoint of a face and edge of respectively. 
However, the vertex solutions do not have an intersection of nodal planes. Instead, the vertex 
solutions have an axis of three- fold symmetry that intersects a vertex in fi. 

The face solutions (type S14) that bifurcate at s = 9 have a bifurcation at s ~ 6.60, as seen in 
Figure [Til The secondary branch (type S51) itself has a bifurcation that creates a tertiary branch 
with type Sgo- Solutions from these new branches are shown in Figure [T6l 

5.3. A degenerate bifurcation with Dg symmetry. Figure [T7] shows the bifurcation diagram 
of the primary branches that bifurcate from the fourth eigenvalue s = \3,i 7 l = -T This is a 
bifurcation with H)q symmetry; the action of H>6 on the critical eigenspace E is shown in Figure [THJ 
Contour plots of the primary branches that bifurcate at s = 11 are shown in Figure [T9l 

The critical eigenspace is the three-dimensional space E = span{V'3,i,i, ^1,3,1 > ^1, 1,3}- The action 
of Tq on E satisfies 

r /r' = {(Rgo, 1)1*, (R 120 , 1)1*, (J 3j -i)r ) D 6 . 
The action of Tq/T' on E = M 3 is isomorphic to the natural action of (M, -R12O) — ^3) on H^ 3 ) where 

r 1 1 
M = 100 . 
.001. 
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Figure 18. The critical eigenspace E at the bifurcation with U)q symmetry of the 
trivial solution at s = 11. The first two figures are similar to those of Figure El 
The arrows show the symmetry of the reduced gradient in E. This is a degenerate 
bifurcation since E is not an irreducible representation space. One of the diagonals of 
the cube in E is an irreducible subspace. The orthogonal subspace, which intersects 
the cube in a hexagon as shown in the third figure, is another irreducible subspace. 

Note that M acts as a reflection. The matrix group (M, i?i2o> — ^3) is called ID^ in the Schonflies 
notation [17] for crystallographic point groups. 

We now describe the symmetries of the solutions shown in Figure [TUl For convenience, we will 
let Tj denote the symmetry of the solution shown with symmetry type Si. 

There is a one-dimensional fixed-point subspace of E for the only symmetry in S2 = {^2}'- 

Fix(r 2 , E) = span{^ 3 ,i,i + "01,3,1 + ^1,1,3}- 

This space is the line through the front and back vertices shown as large dots in the first two cubes 
of Figure [T8l The EBL guarantees that there is a solution with this symmetry; one such branch 
is shown in Figure [T71 The other, negative, branch is not shown. Figure [191 shows a contour map 
of this solution with T2 = symmetry. The solution has one sign on the shaved cube, as shown, 
but the sign is opposite at the center of the cube. The nodal surface has cubic symmetry and is 
diffeomorphic to a sphere. 

Figure [191 also shows one solution with symmetry Tig E S\%. The one-dimensional fixed-point 
subspace of E for this symmetry is 

Fix(ri8,£0 = span{V>3,i,i -"0l,3,i}- 

This fixed-point subspace is the line through the midpoints of two opposite edges, depicted as the 
thickest hexagon diagonal in the third cube of Figure [181 The two other diagonals of the hexagon 
are conjugate fixed-point subspaces. The 180° rotation about each of these diagonals is a symmetry 
of E. There is one conjugacy class of branches that bifurcates with symmetry type Sis at s = 11. 
There are 6 branches in this conjugacy class, one of which is shown in Figure flTl 

Figure fl~9l shows three solutions with symmetry Tn € Si 7. There can be more than one conjugacy 
class of branches because the fixed-point subspace of the Tn action on E is two-dimensional: 

Fix(T 17 ,E) = span{V>3,i,i + "01,3,1, "01,1,3}- 

The intersection of this plane with a cube in E is indicated by dotted lines in Figure [THJ This 
fixed-point subspace includes the one-dimensional intersection of the AIS A\ t 1,3 with E, 

span-^1,1,3} C Fix(T 17 ,E), 

so there is a bifurcating solution branch, which is approximately a multiple of "01,1,3, in .4.1,1,3 ( see 
Equation (jlip ). It is clear from the contour map that the fourth branch from the top has solutions 
that are in .4.1,1,3. Note that this branch undergoes an anomaly-breaking bifurcation at s ~ 3.417, 
with a daughter branch that has the same Tn symmetry. 
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Figure 19. Contour plots of solutions that bifurcate at s = 11. The order of the 
solutions is the same as that in Figure [T71 from top to bottom. The second solution 
branch intersects the fourth solution branch at an anomaly-breaking bifurcation. As 
described in Figure [T71 both have symmetry type S17. The fourth solution's contour 
plot shows a function which is the negative of the continuation of the second solution, 
due to the ink-saving heuristic that replaced u by — u (see Section [3~i|) . 
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Figure 20. Partial bifurcation diagram for PDE ([TJ on the cube, showing the 
primary branch that bifurcates at s = \2,2,2 = 12 and the two bifurcations with 
symmetry that the primary branch undergoes at s ~ 3.697 and s ~ 8.107. There 
is no bifurcation where the branches appear to cross but no circle is drawn. The 
small numerals indicate the MI of the trivial branch, the primary branch, and the 
solutions emanating from the bifurcation with symmetry on the left. Contour 
plots of these solutions are shown in Figure [2TJ 

5.4. A bifurcation with symmetry. Our C++ program can analyze the bifurcations of 
nontrivial solutions, and follow all of the daughter branches of most bifurcations. For example, the 
primary branch that bifurcates at s = 12 undergoes three bifurcations in the interval < s < 12, 
as shown in Figure [20l Two of these bifurcations, at s ~ 3.687 and s ~ 8.107, are bifurcations 
with symmetry. The third bifurcation, at s ~ 8.547, is a bifurcation with B>q symmetry to be 
discussed later. 

We focus on the bifurcation with symmetry at s ~ 3.687. The symmetry of the mother 
solution u* is 

Tl = ((#90, —1), (#120, 1), (#180, — 1), (—h, — 1)}, 
that is, u* G Fix(Ti). For example, a rotation by 90° about the z-axis, coupled with a sign change, 
leaves u* unchanged. The solution on the mother branch shown in Figure I2T1 looks very much like 
^2,2,2- The critical eigenspace E has the ordered basis 

(^2,1,1, ^1,2,1,^1,1,2), 

where ipijk * s a function with the same symmetry as that of ipij k- We cannot find u* or the critical 
eigenfunctions exactly with pencil and paper, but we do know the symmetry exactly, and our C++ 
program is able to use this information. 

The action of T\ on E satisfies = ((—I3, —1)) and the action of Ti/T^ on E is isomorphic to 
the natural action of 

= (— -R90, -R120, — -R180) 

on the coordinate space [E] = M 3 . 

The symmetry of the reduced vector field g on E for this bifurcation with symmetry is shown 
in Figure [3) The daughter solutions at this bifurcation can be classified as face solutions or vertex 
solutions. Each one-dimensional fixed-point subspace of the action on E is conjugate to one of 
these two: 

face [Fix(r 48 , E)} = {(0, 0, a) \ a G M}, vertex [Fix(r 22 , E)} = {(a, 0, 0) | a G M}. 
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M: Si 



V2'. 5*22 



F: S 



18 








Figure 2 1 . Contour plot of the mother solution M shown in Figure [20] along with 
two vertex solutions and a face solution born at a bifurcation with symmetry. 
The names vertex and face indicate the position in the critical eigenspace, shown in 
the middle image of Figure El 



5*98 





Figure 22. Contour plot for a solution with trivial symmetry to PDE ([I]) on the 
cube, at s = 0. This solution is a descendent of the mother branch with symmetry 
type Si shown in Figures [20] and [2T1 The solution shown is found by following the 

sequence of bifurcations dq — > 01 — > 037 — > 091 — > Sqs ■ 
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FIGURE 23. The lattice of isotropy subgroups for the six-dimensional critical 
eigenspace E of the trivial solution at s = 14. For clarity, we display a 3-element 
partition of the edge set of the Hasse diagram of the lattice. The number at the left 
indicates the dimension dim (Fix (r, E)) of the fixed-point subspace for any V € Si 
at that height in the diagram. 
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Figure 24. Symmetry, multiplicity and MI at s = s~ = 14— e of the trivial solution 
and the bifurcating primary solutions at s = 14. The solutions in a given column 
have symmetry type S{. The second row shows the size of the group orbits, while the 
third row gives the MI of the solutions in each group orbit. The only local solution 
at s = s + = 14 + e is the trivial solution, with MI = 17. Using these MI values, we 
can make an index theory computation (|15p to verify that the results are consistent 
with having obtained all solutions. 

Note that I3, the inversion through the origin, is not in T<f. In particular, two antipodal vertex 
solutions are not conjugate, and there is a transcritical branch of vertex solutions, as seen in 
Figure [20] leading to the vertex solutions V\ and V2 seen in Figure EH Note that V\ has two white 
regions and two black regions on the surface of the cube, whereas V2 has one white region and one 
black region. 

Figure [20] shows that the vertex solution V2 is a daughter of both of the bifurcations with 
symmetry on the primary branch that bifurcates at s = 12. The third bifurcation on that branch, 
at s ~ 8.547, is a generic bifurcation with Oq symmetry. The MI of the mother branch changes 
from 11 to 9 as s decreases through that bifurcation. Unlike the degenerate bifurcation with H)q 
symmetry that occurs at u = 0, s = 11, generic bifurcations with Bg symmetry are well-known. 
Hence, we do not give the details of the bifurcation at s ~ 8.547, except to mention that one of the 
two branches created at this bifurcation has a grand-daughter with trivial symmetry, depicted in 
Figure ED 

5.5. A six-dimensional critical eigenspace. Figures [231 EH and [25] concern the bifurcation of 
the trivial solution at s = Ai 2,3 = l 2 + 2 2 + 3 2 = 14. The six-dimensional critical eigenspace is 

E = span{^i )2) 3, "01,3,2, ^2,1,3, "02,3,1, ^3,1,2, ^3,2,l}- 

The action of Tq on E satisfies T' = ((— 13, — 1)) and To/Tp = O^. The action of Tq/T' on E is 
isomorphic to the natural action of 

(i? 90 (-Rsn),Rua Rno, Riso © (--Riso), (~h) © (-I3)) 
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Figure 25. Contour plots for one solution on each of the 19 non-conjugate primary 
branches bifurcating at s = 14. The solutions are listed with increasing MI within 
each symmetry type. The solutions are shown at s = 11, except for 5*67 with MI 
13 and S93, which are shown at s = 13.95 and s = 13.69, respectively. We do this 
because these branches end at s « 13.91 and s ps 13.38, respectively. 
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on the coordinate space [E] = R 6 with respect to the ordered basis 



(V>213 + ^231, ^321 + ^123, ^132 + V>312, V>213 ~ ^231,^321 ~ V>123, V>132 ~ ^312)- 

Writing the action in block diagonal form, one sees that the eigenspace E is the direct sum of 
two irreducible spaces. The trivial subspace of E has isotropy So, and S93 is the minimal isotropy 
subgroup. Thus, [Fix(To, -E)] = {0} C R 6 and [Fix(rg3, E)] = R 6 . The remaining symmetries in 
Figures [23] and [25] in our chosen coordinate space satisfy the following: 



[Fix(r 52 , E)\ = {(0, a, 0, 0, b, 0) | a, b G R}, [Fix(r 54 , E)] = {(a, a, a, b, b,b) \ a,b G R}, 
[Fix(r 44 , E)) = {{a, -a, 0, b, b, 0) | a, b G R}, [Fix(r 78 , E)} = {(a, -a, b, c, -c, 0) | a, b, c G R}, 
[Fix(r 67 , E)} = {(a, b, 0, c, d, 0) \a,b,c,de R}, [Fix(r 79 , £)] = {(a, a, 0, 6, b, c)\a,b,c€ R}. 

Figure [23] describes the lattice of isotropy subgroups of E. Each arrow Si — > Sj indicates that 
some isotropy subgroup in Sj is a subgroup of some isotropy subgroup in Si. The arrows generate 
a partial ordering of the symmetry types. Note that the lattice of isotropy subgroups is different 
from the bifurcation digraph, as explained in |13j . 

To simplify the visual representation, the lattice of symmetry types for the action of Oh on two 
irreducible spaces whose direct sum is E are shown on the top row of Figure [23] Note that the 
middle column is the same in each of the top row sub-lattices. As a result of the presence of S44 
in both sub-lattices, dim(^ n Fix(r 44 )) = 2, whereas dim(^ n Fix^)) = 1 for i G {11, 12, 22, 23}. 
Within these one-dimensional spaces there is a pitchfork bifurcation to an EBL branch, but the 
bifurcation to solutions with symmetry type S* 44 is more complicated. 

It is remarkable that there is at least one solution branch bifurcating at s = 14 with each of the 
symmetry types shown in Figure [23] There is even a solution with symmetry type £93, the lowest 
symmetry present in E. The conjugacy class of this branch has a total of 48 branches. Figure [25] 
shows one solution in each of the nonconjugate primary branches that bifurcate at the multiplicity 
six eigenvalue s = 14. Since each of the solutions in this figure is odd about the center of the cube, 
that is u(x, y, z) = —u(ir — x, ir — y, ir — z), we only show the front view of the contour plot. 

The solution in Figure [25] with symmetry type S52 strongly resembles the eigenfunction 
^i,2,3(x,y, z) = sin(x) sin(2y) sin(3z). An analysis of the "hidden symmetries" in this problem 
|10j would explain why the solution with symmetry type S52 bifurcates, but it would not explain 
all of the solutions in Figure [25l In the space of triply periodic functions on R 3 , there is a 12- 
dimensional irreducible space spanned by rotations of ^1,2,3 and the similar functions with cosines 
in place of sines. 

Let X~ and X + be the set of solutions for s = s~ = 14 — e and s = s + = 14 + e, respectively, that 
are on branches bifurcating from (0, 14) G H x R, together with solutions on the mother branch, for 
a sufficiently small positive e. The set X~ contains 345 solutions falling into 20 group orbits with 
non-trivial representatives shown in Figure [25l Figure [24] shows multiplicity and MI information 
for X~ . Since all the bifurcating branches curve to the left, X + contains only the trivial solution 
with MI(0, s + ) = 17. Thus, one can verify that the Poincare-Hopf Index Theorem of [2] is satisfied 
since 



[Fix(r 12 ,£)] 
[Fix(r 22 ,£)] 



{(0,0, a, 0,0,0) I a G R}, 
{(a, a, a, 0,0,0) | a G R}, 



[Fix(r n ,£)] 
[Fix(r 23 ,£)] 



{(0,0,0, 0,0, a) I a G R}, 
{(0,0,0, a,a,a)\a£ R}, 




This is consistent with our belief that the list of solutions in Figure [Ml is comprehensive. 



6. Conclusion 



In this article we have extended the methods from |12j and [13] . In the first paper, the symmetry 
group was relatively small and a medium-sized grid was used to produce a reasonable portion of the 
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bifurcation diagram and a selection of contour plots for a two-dimensional semilinear elliptic PDE. 
In the second article, we completely automated the symmetry analysis for investigating the rich 
symmetries of solutions to partial difference equations (PdE) for many interesting low-order graphs. 
In the current article, we have shown how to extend these ideas to a three-dimensional problem with 
a large symmetry group, namely the cube. The large grid and many calculations required the use of 
a parallel programming environment. We developed and employed our own library, MPQueue, to 
implement our branch following and branch switching algorithms using self-submitting parallel job 
queues and MPI. The new results we have presented here use the symmetry analysis from a low- 
order graph with the same symmetry group to generate the corresponding symmetry information 
for functions discretized over the large-sized grid used in the PDE code. We could not have done 
this without GAP; for the cube there are 99 symmetry types (and 323 symmetries or isotropy 
subgroups) with 482 arrows between symmetry types. This symmetry information is essential to 
the numerical results in several key ways. It allows for the efficient construction of block-diagonal 
Hessians, reducing the number of costly integrations required at each Newton step in the presence of 
symmetry. It allows us to search for only a single representative of each novel solution type, rather 
than wasting computations on finding many equivalent copies. By reducing the dimension of a 
search space, symmetry information increases our chance of finding all expected solutions of a given 
symmetry type at each new bifurcation. The entire suite of programs and new methods for efficiently 
implementing our algorithms has allowed us to observe interesting bifurcation symmetries for PDE 
that we have not previously seen published. Our procedure demonstrates a robustness for handling 
degenerate bifurcations, AIS, and high-dimensional/reducible critical eigenspaces. The new contour 
plots required a number of ideas for efficiently and effectively conveying the necessary information 
graphically. The size of the problem makes it impossible to present a visual representation of the 
bifurcation digraph on a single page. We have constructed a companion website for navigating the 
digraph, and give examples here to aid the reader in understanding the digraph and how to use it 
to interpret our numerical results. 
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